Namespace Numerical
  ''' <summary>
  ''' Class of helper functions for simple algebra operations on 1 and 2 dimensional double arrays
  ''' Although speed is not essential, we try to avoid recreating and reallocating output arrays 
  ''' on every call as this could slow things down a lot. This means that usually the output arrays MUST
  ''' be allocated and passed to the functions.
  ''' </summary>
  ''' <remarks>One dimensional arrays are called vectors, 2-dimensional arrays are called matrices.
  ''' </remarks>
  Public Class Algebra
    Enum Metric As Integer
      AbsoluteValue = 1
      Euclidean = 2
      Max = 16
    End Enum

    'Sealed class!
    Private Sub New()

    End Sub

    ''' <summary>
    ''' Return a vector with ones
    ''' </summary>
    ''' <param name="iNrCols"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Ones(ByVal iNrCols As Integer) As Double()
      Dim sM(iNrCols - 1) As Double
      Fill(sM, 1)
      Return sM
    End Function

    ''' <summary>
    ''' Return a matrix with ones
    ''' </summary>
    ''' <param name="iNrCols"></param>
    ''' <param name="iNrRows"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Ones(ByVal iNrCols As Integer, ByVal iNrRows As Integer) As Double(,)
      Dim sM(iNrRows - 1, iNrCols - 1) As Double
      Fill(sM, 1)
      Return sM
    End Function

    ''' <summary>
    ''' Return a vector with zeros
    ''' </summary>
    ''' <param name="iNrCols"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Zeros(ByVal iNrCols As Integer) As Double()
      Dim sM(iNrCols - 1) As Double
      Fill(sM, 0)
      Return sM
    End Function

    ''' <summary>
    ''' Return a matrix with ones
    ''' </summary>
    ''' <param name="iNrCols"></param>
    ''' <param name="iNrRows"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Zeros(ByVal iNrCols As Integer, ByVal iNrRows As Integer) As Double(,)
      Dim sM(iNrRows - 1, iNrCols - 1) As Double
      Fill(sM, 0)
      Return sM
    End Function


    ''' <summary>
    ''' Check if each component of a vector is greater than the components of another vector
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <returns></returns>
    ''' <remarks>This method is fastes if the expected result is false.</remarks>
    Public Shared Function IsLargerThan(ByVal v1 As Double(), ByVal v2 As Double()) As Boolean
      For i As Integer = 0 To v1.GetUpperBound(0)
        If v1(i) <= v2(i) Then Return False
      Next
      Return True
    End Function

    ''' <summary>
    ''' Check if each component of a vector is greater than or equal to the components of another vector
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <returns></returns>
    ''' <remarks>This method is fastes if the expected result is false.</remarks>
    Public Shared Function IsLargerThanOrEqual(ByVal v1 As Double(), ByVal v2 As Double()) As Boolean
      For i As Integer = 0 To v1.GetUpperBound(0)
        If v1(i) < v2(i) Then Return False
      Next
      Return True
    End Function


    ''' <summary>
    ''' Check if each component of a vector is smaller than the components of another vector
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <returns></returns>
    ''' <remarks>This method is fastes if the expected result is false.</remarks>
    Public Shared Function IsSmallerThan(ByVal v1 As Double(), ByVal v2 As Double()) As Boolean
      For i As Integer = 0 To v1.GetUpperBound(0)
        If v1(i) >= v2(i) Then Return False
      Next
      Return True
    End Function

    ''' <summary>
    ''' Check if each component of a vector is smaller than the components of another vector
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <returns></returns>
    ''' <remarks>This method is fastes if the expected result is false.</remarks>
    Public Shared Function IsSmallerThanOrEqual(ByVal v1 As Double(), ByVal v2 As Double()) As Boolean
      For i As Integer = 0 To v1.GetUpperBound(0)
        If v1(i) > v2(i) Then Return False
      Next
      Return True
    End Function

    ''' <summary>
    ''' Add 2 vectors
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Add(ByVal sV1() As Double, ByVal sV2() As Double) As Double()
      Dim result(sV1.GetUpperBound(0)) As Double
      Add(sV1, sV2, result)
      Return result
    End Function


    ''' <summary>
    ''' Add 2 vectors
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <param name="sR"></param>
    ''' <remarks>
    ''' </remarks>

    Public Overloads Shared Sub Add(ByVal sV1() As Double, ByVal sV2() As Double, ByVal sR() As Double)
      Dim i, iHiCol As Integer
      iHiCol = sV1.GetUpperBound(0)
      For i = 0 To iHiCol
        sR(i) = sV1(i) + sV2(i)
      Next
    End Sub

    ''' <summary>
    ''' Add 2 matrices together
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Add(ByVal sV1(,) As Double, ByVal sV2(,) As Double) As Double(,)
      Dim result(sV1.GetUpperBound(0), sV1.GetUpperBound(1)) As Double
      Add(sV1, sV2, result)
      Return result
    End Function

    ''' <summary>
    ''' Add 2 matrices
    ''' </summary>
    ''' <param name="sM1"></param>
    ''' <param name="sM2"></param>
    ''' <param name="sMR"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Add(ByVal sM1(,) As Double, ByVal sM2(,) As Double, ByVal sMR(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM1, iHiRow, iHiCol)
      For j = 0 To iHiCol
        For i = 0 To iHiRow
          sMR(i, j) = sM1(i, j) + sM2(i, j)
        Next i
      Next j
    End Sub

    ''' <summary>
    ''' Subtract 2 vectors from each other
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Subtract(ByVal sV1() As Double, ByVal sV2() As Double) As Double()
      Dim result(sV1.GetUpperBound(0)) As Double
      Subtract(sV1, sV2, result)
      Return result
    End Function


    ''' <summary>
    ''' Subtract 2 vectors from each other
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <param name="sR"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Subtract(ByVal sV1() As Double, ByVal sV2() As Double, ByVal sR() As Double)
      Dim i As Integer, iHiCol As Integer
      iHiCol = sV1.GetUpperBound(0)
      For i = 0 To iHiCol
        sR(i) = sV1(i) - sV2(i)
      Next
    End Sub

    ''' <summary>
    '''  Subtract 2 matrices from each other
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Subtract(ByVal sV1(,) As Double, ByVal sV2(,) As Double) As Double(,)
      Dim result(sV1.GetUpperBound(0), sV1.GetUpperBound(1)) As Double
      Subtract(sV1, sV2, result)
      Return result
    End Function


    ''' <summary>
    ''' Subtract 2 matrices from each other
    ''' </summary>
    ''' <param name="sM1"></param>
    ''' <param name="sM2"></param>
    ''' <param name="sMR"></param>
    ''' <remarks>
    ''' </remarks>
    ''' <history>
    ''' 	[yvdh]	25-Oct-05	Created
    ''' </history>

    Public Overloads Shared Sub Subtract(ByVal sM1(,) As Double, ByVal sM2(,) As Double, ByVal sMR(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM1, iHiRow, iHiCol)
      For j = 0 To iHiCol
        For i = 0 To iHiRow
          sMR(i, j) = sM1(i, j) - sM2(i, j)
        Next i
      Next j
    End Sub


    ''' <summary>
    ''' Compute the Euclidean norm of a vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    ''' <history>
    ''' 	[yvdh]	25-Oct-05	Created
    ''' </history>

    Public Overloads Shared Function Norm(ByVal sV1() As Double) As Double
      Return Norm(sV1, Metric.Euclidean)
    End Function


    ''' <summary>
    ''' Compute the Euclidean norm of a 2D vector
    ''' </summary>
    ''' <param name="sX"></param>
    ''' <param name="sY"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    ''' <history>
    ''' 	[yvdh]	25-Oct-05	Created
    ''' </history>

    Public Overloads Shared Function Norm(ByVal sX As Double, ByVal sY As Double) As Double
      Return Math.Sqrt(sX * sX + sY * sY)
    End Function


    ''' <summary>
    ''' Compute the Euclidean norm of a 3D vector
    ''' </summary>
    ''' <param name="sX"></param>
    ''' <param name="sY"></param>
    ''' <param name="sZ"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Norm(ByVal sX As Double, ByVal sY As Double, ByVal sZ As Double) As Double
      Return Math.Sqrt(sX * sX + sY * sY + sZ * sZ)
    End Function


    ''' <summary>
    ''' Compute the N-norm of a vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="iOrder"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Norm(ByVal sV1() As Double, ByVal iOrder As Metric) As Double
      Dim i As Integer, sNorm As Double = 0.0, iHiCol As Integer
      iHiCol = sV1.GetUpperBound(0)
      Select Case iOrder
        Case Metric.AbsoluteValue
          For i = 0 To iHiCol
            sNorm += Math.Abs(sV1(i))
          Next
        Case Metric.Euclidean
          For i = 0 To iHiCol
            sNorm += sV1(i) ^ 2
          Next
          sNorm = sNorm ^ 0.5
        Case Metric.Max
          sNorm = 0
          For i = 0 To iHiCol
            Dim sTemp As Double = Math.Abs(sV1(i))
            If sTemp > sNorm Then sNorm = sTemp
          Next
      End Select
      Return sNorm
    End Function

    ''' <summary>
    ''' Return the Euclidean (N = 2) distance between 2 vectors
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Distance(ByVal v1() As Double, ByVal v2() As Double) As Double
      Return Norm(Subtract(v1, v2))
    End Function

    ''' <summary>
    ''' Return the N-distance between 2 vectors
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <param name="order"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Distance(ByVal v1() As Double, ByVal v2() As Double, ByVal order As Metric) As Double
      Return Norm(Subtract(v1, v2), order)
    End Function

    ''' <summary>
    ''' Compute the columwise mean of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Mean(ByVal sM(,) As Double) As Double()
      Dim sV(sM.GetUpperBound(1)) As Double
      Mean(sM, sV)
      Return sV
    End Function


    ''' <summary>
    ''' Return the mean of a set of points not stored in a matrix
    ''' </summary>
    ''' <param name="datalist"></param>
    ''' <param name="sV"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Mean(ByVal datalist As Generic.IList(Of Double()), ByVal sV() As Double)
      If datalist.Count > 0 Then
        sV.Initialize()
        Dim nrDims As Integer = datalist(0).GetUpperBound(0) + 1
        For Each pt As Double() In datalist
          For i As Integer = 0 To nrDims - 1
            sV(i) += pt(i)
          Next
        Next
        For i As Integer = 0 To nrDims - 1
          sV(i) /= datalist.Count
        Next
      End If
    End Sub


    ''' <summary>
    ''' Return the mean of a set of points not stored in a matrix
    ''' </summary>
    ''' <param name="datalist"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Mean(ByVal datalist As Generic.IList(Of Double())) As Double()
      Dim nrDims As Integer = datalist(0).GetUpperBound(0) + 1
      Dim sV(nrDims - 1) As Double
      Mean(datalist, sV)
      Return sV
    End Function

    ''' <summary>
    ''' Compute the columwise mean of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Mean(ByVal sM(,) As Double, ByVal sV() As Double)
      Sum(sM, sV)
      Multiply(sV, 1.0 / (sM.GetUpperBound(0) + 1), sV)
    End Sub

    ''' <summary>
    ''' Compute the mean of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Mean(ByVal sV() As Double) As Double
      Dim sMean As Double
      sMean = Sum(sV)
      sMean /= sV.GetLength(0)
      Return sMean
    End Function

    ''' <summary>
    ''' Compute the columwise median of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Median(ByVal sM(,) As Double) As Double()
      Dim sV(sM.GetUpperBound(0)) As Double
      Median(sM, sV)
      Return sV
    End Function

    ''' <summary>
    ''' Return the median of each column of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Median(ByVal sM(,) As Double, ByVal sV() As Double)
      Dim i As Integer
      Dim rows As Integer = sM.GetUpperBound(0)
      For i = 0 To rows
        Dim column(rows - 1) As Double
        GetMatrixColumn(sM, column, i)
        sV(i) = Median(column)
      Next
    End Sub

    ''' <summary>
    ''' Return the median of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Median(ByVal sV() As Double) As Double
      'Sort the array
      Array.Sort(sV)
      If sV.Length Mod 2 > 0 Then
        'Odd number of rows
        Return sV(sV.Length \ 2)
      Else
        'Even number of rows
        Return (sV(sV.Length \ 2 - 1) + sV(sV.Length \ 2)) / 2.0
      End If
    End Function

    ''' <summary>
    ''' Compute the columwise percentile of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="p"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Percentile(ByVal sM(,) As Double, ByVal p As Double) As Double()
      Dim sV(sM.GetUpperBound(0)) As Double
      Percentile(sM, p, sV)
      Return sV
    End Function


    ''' <summary>
    ''' Return the Percentile of each column of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="p"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Percentile(ByVal sM(,) As Double, ByVal p As Double, ByVal sV() As Double)
      Dim i As Integer
      Dim rows As Integer = sM.GetUpperBound(0)
      For i = 0 To rows
        Dim column(rows - 1) As Double
        GetMatrixColumn(sM, column, i)
        sV(i) = Percentile(column, p)
      Next
    End Sub


    ''' <summary>
    ''' Return the Percentile of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="p"></param>
    ''' <returns></returns>
    ''' <remarks>Unlike the median, this performs a simple rounding to compute the index of the given percentile.
    ''' </remarks>
    Public Overloads Shared Function Percentile(ByVal sV() As Double, ByVal p As Double) As Double
      'Sort the array
      Array.Sort(sV)
      Dim index As Integer = CInt(Math.Round(sV.Length * p / 100))
      Return sV(index)
    End Function

    ''' <summary>
    ''' Compute the columnwise variances of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <remarks>We use N-1 as normalisation, with N the number of rows in the matrix
    ''' </remarks>
    ''' <history>
    ''' 	[yvdh]	25-Oct-05	Created
    ''' </history>

    Public Overloads Shared Sub Variance(ByVal sM(,) As Double, ByVal sV() As Double)
      Dim iNrCols As Integer = sV.GetUpperBound(0) + 1
      Dim sMean(iNrCols - 1) As Double, sMeanMatrix(sM.GetUpperBound(1), iNrCols - 1) As Double

      Mean(sM, sV)
      Product(Ones(sM.GetUpperBound(1)), sMean, sMeanMatrix)
      Subtract(sM, sMeanMatrix, sM)
      Sum(sM, sV)
      Multiply(sV, sV, sV)
      Multiply(sV, 1.0 / (iNrCols - 1), sV)
    End Sub


    ''' <summary>
    ''' Compute the variance of a vector.
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks>We use N-1 as normalisation, with N the number of elements in the vector
    ''' </remarks>
    Public Overloads Shared Function Variance(ByVal sV() As Double) As Double
      Dim sVar As Double
      Dim sMean(sV.GetUpperBound(0)) As Double
      Fill(sMean, Mean(sV))
      Subtract(sV, sMean, sV)
      Multiply(sV, sV, sV)
      sVar = Sum(sV)
      sVar /= sV.GetUpperBound(0)
      Return sVar
    End Function


    ''' <summary>
    ''' Compute the variance-covariance matrix of a data set in a matrix, with points in rows and components in the columns
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="mean"></param>
    ''' <param name="var"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Covariance(ByVal sV(,) As Double, ByVal mean() As Double, ByVal var(,) As Double)
      Dim nrPts As Integer = sV.GetUpperBound(0) + 1
      Dim nrDims As Integer = sV.GetUpperBound(1) + 1
      Dim temp As Double
      For i As Integer = 0 To nrDims - 1
        For j As Integer = i To nrDims - 1
          temp = 0
          For k As Integer = 0 To nrPts - 1
            temp += (sV(k, i) - mean(i)) * (sV(k, j) - mean(j))
          Next
          temp /= (nrPts - 1)
          var(i, j) = temp
          If i <> j Then var(j, i) = temp
        Next
      Next
    End Sub

    ''' <summary>
    ''' Compute the variance-covariance matrix of a data set in a matrix, with points in rows and components in the columns
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Function Covariance(ByVal sV As Generic.IList(Of Double())) As Double(,)
      Dim nrPts As Integer = sV.Count
      Dim nrDims As Integer = sV(0).GetUpperBound(0) + 1
      Dim var(nrDims - 1, nrDims - 1) As Double
      Dim avg() As Double = Mean(sV)
      Dim temp As Double
      For i As Integer = 0 To nrDims - 1
        For j As Integer = i To nrDims - 1
          temp = 0
          For Each pt As Double() In sV
            temp += (pt(i) - avg(i)) * (pt(j) - avg(j))
          Next
          temp /= (nrPts - 1)
          var(i, j) = temp
          If i <> j Then var(j, i) = temp
        Next
      Next
      Return var
    End Function

    ''' <summary>
    ''' Compute the columwise sum of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Sum(ByVal sM(,) As Double) As Double()
      Dim sV(sM.GetUpperBound(0)) As Double
      Sum(sM, sV)
      Return sV
    End Function


    ''' <summary>
    ''' Compute columnwise sum of matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Sum(ByVal sM(,) As Double, ByVal sV() As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM, iHiRow, iHiCol)

      For j = 0 To iHiCol
        sV(j) = 0.0
        For i = 0 To iHiRow
          sV(j) += sM(i, j)
        Next i
      Next j
    End Sub


    ''' <summary>
    ''' Compute sum of elements of vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Sum(ByVal sV() As Double) As Double

      Dim sSum As Double = 0
      Dim i, iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)

      For i = 0 To iHiCol
        sSum = sSum + sV(i)
      Next i
      Return sSum
    End Function


    ''' <summary>
    ''' Find the maximum value of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="iPos">Contains the index of the highest value</param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Max(ByVal sV() As Double, ByRef iPos As Integer) As Double
      Dim i As Integer, sMax As Double = Double.MinValue
      Dim iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)

      For i = 0 To iHiCol
        If sV(i) > sMax Then
          iPos = i
          sMax = sV(i)
        End If
      Next i
      Return sMax
    End Function


    ''' <summary>
    ''' Find the maximum value of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Max(ByVal sV() As Double) As Double
      Dim iPos As Integer
      Return Max(sV, iPos)
    End Function


    ''' <summary>
    ''' Find the maximum value of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Max(ByVal sM(,) As Double) As Double
      Dim i, j As Integer
      Return Max(sM, i, j)
    End Function


    ''' <summary>
    ''' Find the maximum value of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iCol">The column of the maximum</param>
    ''' <param name="iRow">The row of the maximum</param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Max(ByVal sM(,) As Double, ByRef iCol As Integer, ByRef iRow As Integer) As Double
      Dim sMAx As Double = Double.MinValue
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM, iHiRow, iHiCol)
      For j = 0 To iHiCol
        For i = 0 To iHiRow
          If sM(i, j) > sMAx Then
            iCol = j
            iRow = i
            sMAx = sM(i, j)
          End If
        Next i
      Next j
      Return sMAx
    End Function


    ''' <summary>
    ''' Find the minimum value of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="iPos">Contains the index of the lowest value</param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Min(ByVal sV() As Double, ByRef iPos As Integer) As Double
      'Find Min of a vector
      Dim i As Integer, sMin As Double = Double.MaxValue
      Dim iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)

      For i = 0 To iHiCol
        If sV(i) < sMin Then
          iPos = i
          sMin = sV(i)
        End If
      Next i
      Return sMin
    End Function


    ''' <summary>
    ''' Find the minimum value of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Min(ByVal sV() As Double) As Double
      Dim iPos As Integer
      Return Min(sV, iPos)
    End Function


    ''' <summary>
    ''' Find the minimum value of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Min(ByVal sM(,) As Double) As Double
      'Find Min of a matrix
      Dim i, j As Integer
      Return Min(sM, i, j)
    End Function


    ''' <summary>
    ''' Find the minimum value of a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iCol">The column of the  minimum</param>
    ''' <param name="iRow">The row of the  minimum</param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Min(ByVal sM(,) As Double, ByRef iCol As Integer, ByRef iRow As Integer) As Double
      'Find Min of a matrix
      Dim sMin As Double = Double.MaxValue
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM, iHiRow, iHiCol)
      For j = 0 To iHiCol
        For i = 0 To iHiRow
          If sM(i, j) < sMin Then
            iCol = j
            iRow = i
            sMin = sM(i, j)
          End If
        Next i
      Next j
      Return sMin
    End Function

    ''' <summary>
    ''' Compute the power of the elements of a vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="exponent"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Pow(ByVal sV1() As Double, _
                                         ByVal exponent As Double) As Double()
      Dim sV2(sV1.GetUpperBound(0)) As Double
      Pow(sV1, exponent, sV2)
      Return sV2
    End Function

    ''' <summary>
    ''' Compute the power of the elements of a vector.
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="exponent"></param>
    ''' <param name="sV2"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Pow(ByVal sV1() As Double, ByVal exponent As Double, ByVal sV2() As Double)
      Dim i As Integer, iHiRow As Integer
      iHiRow = UBound(sV1)
      For i = 0 To iHiRow
        sV2(i) = Math.Pow(sV1(i), exponent)
      Next i
    End Sub

    ''' <summary>
    ''' Compute the power of the elements of a matrix
    ''' </summary>
    ''' <param name="m1 "></param>
    ''' <param name="exponent"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Pow(ByVal m1(,) As Double, _
                                         ByVal exponent As Double) As Double(,)
      Dim m2(m1.GetUpperBound(0), m1.GetUpperBound(1)) As Double
      Pow(m1, exponent, m2)
      Return m2
    End Function

    ''' <summary>
    ''' Compute the power of the elements of a matrix.
    ''' </summary>
    ''' <param name="m1 "></param>
    ''' <param name="exponent"></param>
    ''' <param name="m2 "></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Pow(ByVal m1(,) As Double, ByVal exponent As Double, ByVal m2(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(m1, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          m2(i, j) = Math.Pow(m1(i, j), exponent)
        Next j
      Next i
    End Sub




    ''' <summary>
    ''' Compute the negation of the elements of a vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Negate(ByVal sV1() As Double) As Double()
      Dim sV2(sV1.GetUpperBound(0)) As Double
      Negate(sV1, sV2)
      Return sV2
    End Function

    ''' <summary>
    ''' Compute the negation of the elements of a vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Negate(ByVal sV1() As Double, ByVal sV2() As Double)
      Dim i As Integer, iHiRow As Integer
      iHiRow = UBound(sV1)
      For i = 0 To iHiRow
        sV2(i) = -sV1(i)
      Next i
    End Sub

    ''' <summary>
    ''' Compute the negation  of the elements of a matrix
    ''' </summary>
    ''' <param name="m1 "></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Negate(ByVal m1(,) As Double) As Double(,)
      Dim m2(m1.GetUpperBound(0), m1.GetUpperBound(1)) As Double
      Negate(m1, m2)
      Return m2
    End Function

    ''' <summary>
    ''' Compute the negation  of the elements of a matrix.
    ''' </summary>
    ''' <param name="m1 "></param>
    ''' <param name="m2 "></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Negate(ByVal m1(,) As Double, ByVal m2(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(m1, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          m2(i, j) = -m1(i, j)
        Next j
      Next i
    End Sub


    ''' <summary>
    ''' Multiply a scalar with an offset. For vectors and matrices this would lead to too many 
    ''' different versions, so use Subtract to have an offset.
    ''' </summary>
    ''' <param name="sX"></param>
    ''' <param name="sOffset"></param>
    ''' <param name="sScale"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Multiply(ByVal sX As Double, ByVal sOffset As Double, ByVal sScale As Double) As Double
      Return (sX - sOffset) * sScale
    End Function

    ''' <summary>
    ''' Multiply elements of a vector with a scalar
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sScale"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Multiply(ByVal sV1() As Double, ByVal sScale As Double) As Double()
      Dim sY(sV1.GetUpperBound(0)) As Double
      Multiply(sV1, sScale, sY)
      Return sY
    End Function


    ''' <summary>
    ''' Multiply elements of a vector with a scalar
    ''' </summary>
    ''' <param name="sScale"></param>
    ''' <param name="sV2"></param>
    ''' <param name="sY"></param>
    ''' <remarks>Output vector maybe the same as the input vector
    ''' </remarks>
    Public Overloads Shared Sub Multiply(ByVal sV2() As Double, _
                                         ByVal sScale As Double, _
                                         ByVal sY() As Double)
      Dim i As Integer, iHiRow As Integer
      iHiRow = UBound(sV2)
      For i = 0 To iHiRow
        sY(i) = sScale * sV2(i)
      Next i
    End Sub

    ''' <summary>
    ''' Multiply the elements of 2 vectors
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Multiply(ByVal sV1() As Double, _
                                              ByVal sV2() As Double) As Double()
      Dim sY(sV1.GetUpperBound(0)) As Double
      Multiply(sV1, sV2, sY)
      Return sY
    End Function


    ''' <summary>
    ''' Multiply the elements of 2 vectors
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <param name="sY">The resulting vector. Maybe the same as one of the input vectors</param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Multiply(ByVal sV1() As Double, _
                                      ByVal sV2() As Double, _
                                      ByVal sY() As Double)
      Dim i As Integer, iHiRow As Integer
      iHiRow = UBound(sV2)
      For i = 0 To iHiRow
        sY(i) = sV1(i) * sV2(i)
      Next i
    End Sub


    ''' <summary>
    ''' Multiply elements of a matrix with a scalar
    ''' </summary>
    ''' <param name="sB"></param>
    ''' <param name="sScale"></param>
    ''' <param name="sY">The resulting matrix. Mabye the same as the input matrix</param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Multiply(ByVal sB(,) As Double, _
                                         ByVal sScale As Double, _
                                         ByVal sY(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sY(i, j) = sScale * sB(i, j)
        Next j
      Next i
    End Sub


    ''' <summary>
    ''' Multiply the elements of a matrix using the elements of another matrix
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sB"></param>
    ''' <param name="sY">The resulting matrix. Mabye the same as one of the input matrices</param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Multiply(ByVal sA(,) As Double, _
                                    ByVal sB(,) As Double, _
                                    ByVal sY(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sY(i, j) = sA(i, j) * sB(i, j)
        Next j
      Next i
    End Sub

    ''' <summary>
    ''' Multiply elements of a matrix using the corresponding elements of a vector, per row or per column.
    ''' </summary>
    ''' <param name="sB"></param>
    ''' <param name="scales"></param>
    ''' <param name="row">Use the scales in the vector per row if true</param>
    ''' <param name="sY"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Multiply(ByVal sB(,) As Double, _
                                      ByVal scales() As Double, _
                                      ByVal row As Boolean, _
                                      ByVal sY(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)

      If row Then
        For i = 0 To iHiRow
          For j = 0 To iHiCol
            sY(i, j) = scales(i) * sB(i, j)
          Next j
        Next i
      Else
        For i = 0 To iHiRow
          For j = 0 To iHiCol
            sY(i, j) = scales(j) * sB(i, j)
          Next j
        Next i
      End If
    End Sub

    ''' <summary>
    ''' Divide elements of a vector by a scalar
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sScale"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Divide(ByVal sV1() As Double, _
                                            ByVal sScale As Double) As Double()
      Dim sY() As Double = Multiply(sV1, sScale)
      Return sY
    End Function


    ''' <summary>
    ''' Divide elements of vector by a scalar
    ''' </summary>
    ''' <param name="sV2"></param>
    ''' <param name="sScale"></param>
    ''' <param name="sY">The resulting vector. Mabye the same vector as the input vector</param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Divide(ByVal sV2() As Double, _
                                       ByVal sScale As Double, _
                                       ByVal sY() As Double)
      Multiply(sV2, 1 / sScale, sY)
    End Sub


    ''' <summary>
    ''' Divide elements of a vector using the elements of another vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Divide(ByVal sV1() As Double, _
                                            ByVal sV2() As Double) As Double()
      Dim sY(sV1.GetUpperBound(0)) As Double
      Divide(sV1, sV2, sY)
      Return sY
    End Function


    ''' <summary>
    ''' Divide elements of a vector using the elements of another vector
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <param name="sY">The resulting vector. Mabye the same vector as one of the input vectors</param>
    ''' <remarks>V2 must not contain zeros
    ''' </remarks>
    Public Overloads Shared Sub Divide(ByVal sV1() As Double, _
                                      ByVal sV2() As Double, _
                                      ByVal sY() As Double)
      Dim i As Integer, iHiRow As Integer
      iHiRow = UBound(sV2)
      For i = 0 To iHiRow
        sY(i) = sV1(i) / sV2(i)
      Next i
    End Sub


    ''' <summary>
    ''' Divide elements of a matrix by a scalar
    ''' </summary>
    ''' <param name="sB"></param>
    ''' <param name="sScale"></param>
    ''' <param name="sY">The resulting matrix. Maybe the same as the input matrix</param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Divide(ByVal sB(,) As Double, _
                                       ByVal sScale As Double, _
                                       ByVal sY(,) As Double)
      Multiply(sB, 1 / sScale, sY)
    End Sub


    ''' <summary>
    ''' Divide elements of a matrix using the corresponding elements of another matrix
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sB"></param>
    ''' <param name="sY">The resulting matrix. Maybe the same as one of the input matrices</param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Divide(ByVal sA(,) As Double, _
                                    ByVal sB(,) As Double, _
                                    ByVal sY(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sY(i, j) = sA(i, j) / sB(i, j)
        Next j
      Next i
    End Sub


    ''' <summary>
    ''' Divide elements of a matrix using the corresponding elements of a vector, per row or per column
    ''' </summary>
    ''' <param name="sB"></param>
    ''' <param name="scales"></param>
    ''' <param name="row">Use the scales in the vector per row if true</param>
    ''' <param name="sY"></param>
    ''' <remarks>No zeros in the scales of course!
    ''' </remarks>
    Public Overloads Shared Sub Divide(ByVal sB(,) As Double, _
                                      ByVal scales() As Double, _
                                      ByVal row As Boolean, _
                                      ByVal sY(,) As Double)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)

      If row Then
        For i = 0 To iHiRow
          For j = 0 To iHiCol
            sY(i, j) = sB(i, j) / scales(i)
          Next j
        Next i
      Else
        For i = 0 To iHiRow
          For j = 0 To iHiCol
            sY(i, j) = sB(i, j) / scales(j)
          Next j
        Next i
      End If
    End Sub


    ''' <summary>
    ''' Compute matrix product
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sB"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Product(ByVal sA(,) As Double, _
                                        ByVal sB(,) As Double) As Double(,)
      Dim sC(sA.GetUpperBound(0), sB.GetUpperBound(1)) As Double
      Product(sA, sB, sC)
      Return sC
    End Function

    ''' <summary>
    ''' Compute matrix product
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sB"></param>
    ''' <param name="sC"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Product(ByVal sA(,) As Double, _
                                        ByVal sB(,) As Double, _
                                        ByVal sC(,) As Double)
      'Compute A * B and store in C. 
      'Raise a fatal run-time error if any errors (no return value)!
      'Comment Oct. 2003
      ' A profiler has demonstrated the hefty price of bounds checking in this routine, 
      ' consequently I have commented it!
      'Dim i, j, k, iAHiRow, iAHiCol As Integer
      'GetBounds(sA, iAHiRow, iAHiCol)
      'Dim iBHiRow, iBHiCol As Integer
      'GetBounds(sB, iBHiRow, iBHiCol)
      'Dim iCHiRow, iCHiCol As Integer
      'GetBounds(sC, iCHiRow, iCHiCol)
      'If (((iAHiCol) <> (iBHiRow)) Or _
      '    ((iAHiRow) <> (iCHiRow)) Or _
      '    ((iBHiCol) <> (iCHiCol))) Then
      '  MsgBox("Algebra.Product: Incompatible matrix dimensions", MsgBoxStyle.OKOnly + MsgBoxStyle.Critical)
      'End If

      Dim i, j, k As Integer, iAHiCol As Integer = sA.GetUpperBound(1)
      Dim iCHiRow As Integer = sC.GetUpperBound(0), iCHiCol As Integer = sC.GetUpperBound(1)

      For i = 0 To iCHiRow
        For j = 0 To iCHiCol
          sC(i, j) = 0.0
          For k = 0 To iAHiCol
            sC(i, j) += sA(i, k) * sB(k, j)
          Next k
        Next j
      Next i
    End Sub


    ''' <summary>
    ''' Round a double vector to an integer vector
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function Round(ByVal v1() As Double) As Integer()
      Dim v2(v1.GetUpperBound(0)) As Integer
      Round(v1, v2)
      Return v2
    End Function

    ''' <summary>
    '''  Round a double vector to an integer vector
    ''' </summary>
    ''' <param name="v1"></param>
    ''' <param name="v2"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Round(ByVal v1() As Double, ByVal v2() As Integer)
      Dim hi As Integer = v1.GetUpperBound(0)
      For i As Integer = 0 To hi
        v2(i) = CInt(Math.Round(v1(i)))
      Next
    End Sub

    ''' <summary>
    ''' Return the product of 2 vectors (dot product)
    ''' </summary>
    ''' <param name="sV1"></param>
    ''' <param name="sV2"></param>
    ''' <returns></returns>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Function Product(ByVal sV1() As Double, ByVal sV2() As Double) As Double
      'Return the scalar product of two vectors.
      Dim i As Integer, iHiRow As Integer, sResult As Double

      iHiRow = sV1.GetUpperBound(0)
      For i = 0 To iHiRow
        sResult = sResult + sV1(i) * sV2(i)
      Next i
      Return sResult
    End Function


    ''' <summary>
    ''' Multiply a vector times a vector (Y = M.Y), by interpreting the vector M as a columnmatrix,
    ''' and X as a rowmatrix. Result is a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sX"></param>
    ''' <param name="sY"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Product(ByVal sM() As Double, _
                                        ByVal sX() As Double, _
                                        ByVal sY(,) As Double)

      Dim sA(0, sM.GetUpperBound(0)) As Double, sB(sX.GetUpperBound(0), 0) As Double

      SetMatrixColumn(sA, sM, 0)
      SetMatrixRow(sB, sX, 0)
      Product(sA, sB, sY)
    End Sub


    ''' <summary>
    ''' Multiply a matrix times a vector (y = M.x), by interpreting the vector X as a columnmatrix.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sX"></param>
    ''' <param name="sY"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Product(ByVal sM(,) As Double, _
                                        ByVal sX() As Double, _
                                        ByVal sY() As Double)
      '
      ' COMMENT Oct. 2003: A profiler has shown a hefty price for the 2 temporary matrices when 
      ' using the easy way of coding, i.e. based on another version of product.
      'Dim iHiRow As Integer = sX.GetUpperBound(0)
      'Dim sB(iHiRow, 0), sC(iHiRow, 0) As Double

      'SetMatrixColumn(sB, sX, 0)
      'Product(sM, sB, sC)
      'GetMatrixColumn(sC, sY, 0)

      ' Faster method without temp matrices (5 x faster!)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM, iHiRow, iHiCol)
      For i = 0 To iHiRow
        sY(i) = 0
        For j = 0 To iHiCol
          sY(i) += sM(i, j) * sX(j)
        Next
      Next

    End Sub


    ''' <summary>
    ''' Multiply a vector with a matrix (y = x.M), by interpreting the vector X as a rowmatrix.
    ''' </summary>
    ''' <param name="sX"></param>
    ''' <param name="sM"></param>
    ''' <param name="sY"></param>
    ''' <remarks>
    ''' </remarks>
    Public Overloads Shared Sub Product(ByVal sX() As Double, _
                                        ByVal sM(,) As Double, _
                                        ByVal sY() As Double)
      '
      'See comment at overloaded function for speed questions

      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM, iHiRow, iHiCol)
      For i = 0 To iHiCol
        sY(i) = 0
        For j = 0 To iHiRow
          sY(i) += sX(j) * sM(j, i)
        Next
      Next

    End Sub

    ''' <summary>
    ''' Multiply each element from the list by the matrix M
    ''' </summary>
    ''' <param name="M"></param>
    ''' <param name="dataList"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Product(ByVal M(,) As Double, ByVal dataList As Generic.IList(Of Double()))
      For Each pt As Double() In dataList
        Product(M, pt, pt)
      Next
    End Sub


    ''' <summary>
    ''' Extract a submatrix.
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sB">The resulting submatrix. Set its dimensions beforehand</param>
    ''' <param name="iRow">The starting row</param>
    ''' <param name="iCol">The starting column</param>
    ''' <remarks>
    ''' </remarks>
    Public Shared Sub GetSubMatrix(ByVal sA(,) As Double, _
                         ByVal sB(,) As Double, _
                         ByVal iRow As Integer, _
                         ByVal iCol As Integer)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sB(i, j) = sA(i + iRow, j + iCol)
        Next j
      Next i
    End Sub

    ''' <summary>
    ''' Set a submatrix.
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sB">The submatrix.</param>
    ''' <param name="iRow">The starting row</param>
    ''' <param name="iCol">The starting column</param>
    ''' <remarks>
    ''' </remarks>
    Public Shared Sub SetSubMatrix(ByVal sA(,) As Double, _
                         ByVal sB(,) As Double, _
                         ByVal iRow As Integer, _
                         ByVal iCol As Integer)
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sB, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sA(i + iRow, j + iCol) = sB(i, j)
        Next j
      Next i
    End Sub

    ''' <summary>
    ''' Get a matrix column.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iCol"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function GetMatrixColumn(ByVal sM(,) As Double, ByVal iCol As Integer) As Double()
      Dim sv(sM.GetUpperBound(0)) As Double
      GetMatrixColumn(sM, sv, iCol, 0)
      Return sv
    End Function


    ''' <summary>
    ''' Get a matrix column.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iCol"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub GetMatrixColumn(ByVal sM(,) As Double, _
                            ByVal sV() As Double, _
                            ByVal iCol As Integer)
      GetMatrixColumn(sM, sV, iCol, 0)
    End Sub

    ''' <summary>
    ''' Get a matrix column starting from a given row.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iCol"></param>
    ''' <param name="iStartRow"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function GetMatrixColumn(ByVal sM(,) As Double, _
                            ByVal iCol As Integer, _
                            ByVal iStartRow As Integer) As Double()
      Dim sv(sM.GetUpperBound(0) - iStartRow) As Double
      GetMatrixColumn(sM, sv, iCol, iStartRow)
      Return sv
    End Function


    ''' <summary>
    ''' Get a matrix column starting from a given row.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iCol"></param>
    ''' <param name="iStartRow"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub GetMatrixColumn(ByVal sM(,) As Double, _
                            ByVal sV() As Double, _
                            ByVal iCol As Integer, _
                            ByVal iStartRow As Integer)
      'Fill vector with matrix col
      Dim i, iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)
      For i = 0 To iHiCol
        sV(i) = sM(i + iStartRow, iCol)
      Next i
    End Sub

    ''' <summary>
    ''' Get a matrix row.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iCol"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function GetMatrixRow(ByVal sM(,) As Double, ByVal iCol As Integer) As Double()
      Dim sv(sM.GetUpperBound(1)) As Double
      GetMatrixRow(sM, sv, iCol, 0)
      Return sv
    End Function

    ''' <summary>
    ''' Get a matrix row.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iRow"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub GetMatrixRow(ByVal sM(,) As Double, _
                          ByVal sV() As Double, _
                          ByVal iRow As Integer)
      GetMatrixRow(sM, sV, iRow, 0)
    End Sub

    ''' <summary>
    ''' Get a matrix row starting from a given column.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iRow "></param>
    ''' <param name="iStartCol "></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function GetMatrixRow(ByVal sM(,) As Double, _
                            ByVal iRow As Integer, _
                            ByVal iStartCol As Integer) As Double()
      Dim sv(sM.GetUpperBound(1) - iStartCol) As Double
      GetMatrixRow(sM, sv, iRow, iStartCol)
      Return sv
    End Function

    ''' <summary>
    ''' Get a matrix row starting from a given column.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iRow"></param>
    ''' <param name="iStartCol"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub GetMatrixRow(ByVal sM(,) As Double, _
                            ByVal sV() As Double, _
                            ByVal iRow As Integer, _
                            ByVal iStartCol As Integer)
      Dim i, iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)
      For i = 0 To iHiCol
        sV(i) = sM(iRow, i + iStartCol)
      Next i
    End Sub

    ''' <summary>
    ''' Fill a matrix with a value
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sVal"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Fill(ByVal sM(,) As Double, ByVal sVal As Double)
      'Fill with given number
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sM, iHiRow, iHiCol)
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sM(i, j) = sVal
        Next
      Next
    End Sub

    ''' <summary>
    ''' Fill a vector with a value
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sVal"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Fill(ByVal sM() As Double, ByVal sVal As Double)
      'Fill with given number
      Dim i, iHiCol As Integer
      iHiCol = sM.GetUpperBound(0)
      For i = 0 To iHiCol
        sM(i) = sVal
      Next
    End Sub

    ''' <summary>
    ''' Set a matrix column using a vector, starting from a given row.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iCol"></param>
    ''' <param name="iStartRow"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixColumn(ByVal sM(,) As Double, _
                              ByVal sV() As Double, _
                              ByVal iCol As Integer, _
                              ByVal iStartRow As Integer)
      'Fill matrix col with vector
      Dim i, iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)
      For i = 0 To iHiCol
        sM(i + iStartRow, iCol) = sV(i)
      Next i
    End Sub

    ''' <summary>
    ''' Set a matrix column using a vector
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iCol"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixColumn(ByVal sM(,) As Double, _
                               ByVal sV() As Double, _
                               ByVal iCol As Integer)
      SetMatrixColumn(sM, sV, iCol, 0)
    End Sub

    ''' <summary>
    ''' Set a matrix column using a scalar, starting from a given row.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iCol"></param>
    ''' <param name="iStartRow"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixColumn(ByVal sM(,) As Double, _
                              ByVal sV As Double, _
                              ByVal iCol As Integer, _
                              ByVal iStartRow As Integer)
      'Fill matrix col with vector
      Dim i, iHiCol As Integer
      iHiCol = sM.GetUpperBound(1)
      For i = 0 To iHiCol
        sM(i + iStartRow, iCol) = sV
      Next i
    End Sub

    ''' <summary>
    ''' Set a matrix column using a scalar.
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iCol"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixColumn(ByVal sM(,) As Double, _
                               ByVal sV As Double, _
                               ByVal iCol As Integer)
      SetMatrixColumn(sM, sV, iCol, 0)
    End Sub

    ''' <summary>
    ''' Set a matrix row using a vector
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iRow"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixRow(ByVal sM(,) As Double, _
                          ByVal sV() As Double, _
                          ByVal iRow As Integer)
      SetMatrixRow(sM, sV, iRow, 0)
    End Sub


    ''' <summary>
    ''' Set a matrix row using a vector, starting from a given column
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iRow"></param>
    ''' <param name="iStartCol"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixRow(ByVal sM(,) As Double, _
                            ByVal sV() As Double, _
                            ByVal iRow As Integer, _
                            ByVal iStartCol As Integer)
      Dim i, iHiCol As Integer
      iHiCol = sV.GetUpperBound(0)
      For i = 0 To iHiCol
        sM(iRow, i + iStartCol) = sV(i)
      Next i
    End Sub

    ''' <summary>
    '''  Set a matrix row using a scalar
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iRow"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixRow(ByVal sM(,) As Double, _
                          ByVal sV As Double, _
                          ByVal iRow As Integer)
      SetMatrixRow(sM, sV, iRow, 0)
    End Sub

    ''' <summary>
    '''  Set a matrix row using a vector, starting from a given column
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <param name="iRow"></param>
    ''' <param name="iStartCol"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub SetMatrixRow(ByVal sM(,) As Double, _
                            ByVal sV As Double, _
                            ByVal iRow As Integer, _
                            ByVal iStartCol As Integer)
      Dim i, iHiCol As Integer
      iHiCol = sM.GetUpperBound(0)
      For i = 0 To iHiCol
        sM(iRow, i + iStartCol) = sV
      Next i
    End Sub

    ''' <summary>
    ''' Create a vector from a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function MatrixToVector(ByVal sM(,) As Double) As Double()
      'Put all elements of a matrix into a vector
      Dim V(sM.GetLength(0) * sM.GetLength(1) - 1) As Double
      MatrixToVector(sM, V)
      Return V
    End Function

    ''' <summary>
    ''' Fill a vector from the elements of a matrix. 
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="sV"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub MatrixToVector(ByVal sM(,) As Double, _
                                     ByVal sV() As Double)
      'Put all elements of a matrix into a vector
      Dim i, j, iHiRow, iHiCol, k As Integer

      GetBounds(sM, iHiRow, iHiCol)
      k = 0
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sV(k) = sM(i, j)
          k += 1
        Next
      Next
    End Sub

    ''' <summary>
    ''' Fill the elements of a matrix from those of a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function VectorToMatrix(ByVal sV() As Double, ByVal rows As Integer, ByVal cols As Integer) As Double(,)
      Dim result(rows - 1, cols - 1) As Double
      VectorToMatrix(sV, result)
      Return result
    End Function

    ''' <summary>
    ''' Fill the elements of a matrix from those of a vector.
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Shared Sub VectorToMatrix(ByVal sV() As Double, _
                                     ByVal sM(,) As Double)
      'Put all elements of a vector into a vector. Use the shape of the matrix
      Dim i, j, iHiRow, iHiCol, k As Integer

      GetBounds(sM, iHiRow, iHiCol)
      k = 0
      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sM(i, j) = sV(k)
          k += 1
        Next
      Next
    End Sub

    ''' <summary>
    ''' Transpose a matrix
    ''' </summary>
    ''' <param name="sA"></param>
    ''' <param name="sAt"></param>
    ''' <remarks></remarks>
    Public Shared Sub Transpose(ByVal sA(,) As Double, ByVal sAt(,) As Double)
      'Transpose matrix A and put result in At. Output has
      'same base as input. Input arguments must be different!
      Dim i, j, iHiRow, iHiCol As Integer
      GetBounds(sA, iHiRow, iHiCol)

      For i = 0 To iHiRow
        For j = 0 To iHiCol
          sAt(j, i) = sA(i, j)
        Next j
      Next i
    End Sub

    ''' <summary>
    ''' Read a tex file with a matrix or vector stored separated by spaces and
    '''newlines. sM will be redimensioned as necessary and must be
    '''a dynamic array. Redimensioning can only affect the last dimension!
    '''When a vector is read in the matrix will be of size n x 1, and can easily 
    '''be converted to a vector
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Shared Function Load(ByVal strFile As String) As Double(,)
      Dim sM(,) As Double
      Dim iNrCols As Integer
      Dim iRowNr As Integer, iColNr As Integer
      Dim sMt(,) As Double, strText As String, strTextItems() As String
      Dim objCulture As New System.Globalization.CultureInfo("") 'Use invariant culture!!!!!

      Debug.WriteLine("Loading numeric data file " & strFile)
      Try
        FileOpen(5, strFile, OpenMode.Input, OpenAccess.Read)
      Catch e As Exception
        MsgBox("Algebra.Load:" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
        Return Nothing
      End Try

      iRowNr = 0
      iNrCols = 0
      Do While Not EOF(5)
        'Read first line to count number of columns
        strText = Trim(LineInput(5))

        If strText.Length > 0 Then
          strText = strText.Replace("  ", " ") 'Make sure no 2 spaces are in the string ...
          strText = strText.Replace("   ", " ") 'Make sure no 3 spaces are in the string ...
          strTextItems = strText.Split()

          'Redimension the array if the nr of cols is known, i.e. after
          'reading the first line.
          If iRowNr = 0 Then
            iNrCols = (strTextItems.GetUpperBound(0) + 1)
            ReDim sMt(iNrCols - 1, 0)
          Else
            ReDim Preserve sMt(iNrCols - 1, iRowNr)
          End If

          'Read values into transposed matrix
          For iColNr = 0 To iNrCols - 1
            'sMt(iColNr, iRowNr) = CSng(strTextItems(iColNr))
            sMt(iColNr, iRowNr) = Double.Parse(strTextItems(iColNr), objCulture)
          Next
          iRowNr += 1
        End If
      Loop

      'close file
      FileClose(5)

      'Transpose matrix to output format
      ReDim sM(iRowNr - 1, iNrCols - 1)
      Transpose(sMt, sM)
      Return sM
    End Function

    ''' <summary>
    ''' Save a matrix to a text file
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                        ByVal sM(,) As Double)
      Save(strFile, sM, 16, 2)
    End Sub

    ''' <summary>
    ''' Save a matrix to a text file, using a given precision
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <param name="iPrecBeforeDec"></param>
    ''' <param name="iPrecAfterDec"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                          ByVal sM(,) As Double, _
                          ByVal iPrecBeforeDec As Integer, _
                          ByVal iPrecAfterDec As Integer)
      'Save a matrix to file.
      Dim objCulture As New System.Globalization.CultureInfo("") 'Use invariant culture!!!!!
      Dim strF As String = ""
      Dim i, j, iHiRow, iHiCol As Integer

      If iPrecAfterDec = -1 Then
        strF = "0."
      Else
        For i = 1 To iPrecAfterDec
          strF = strF & "0"
        Next
        strF = strF & "."
      End If

      For i = 1 To iPrecBeforeDec
        strF = strF & "#"
      Next

      If System.IO.File.Exists(strFile) Then System.IO.File.Delete(strFile)

      Try
        FileOpen(5, strFile, OpenMode.Output, OpenAccess.Write)
        GetBounds(sM, iHiRow, iHiCol)
        For i = 0 To iHiRow
          For j = 0 To iHiCol - 1
            'Print(5, Format(sM(i, j), strF), SPC(1))
            Print(5, sM(i, j).ToString(strF, objCulture), SPC(1))
          Next j
          'PrintLine(5, SPC(1), Format(sM(i, iHiCol), strF))
          PrintLine(5, SPC(1), sM(i, iHiCol).ToString(strF, objCulture))
        Next i
        FileClose(5)
      Catch e As Exception
        MsgBox("Algebra.Save (file = " & strFile & "):" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
      End Try
    End Sub

    ''' <summary>
    ''' Save a vector to a file
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                                     ByVal sM() As Double)
      Save(strFile, sM, 16, 2)
    End Sub

    ''' <summary>
    ''' Save a vector to a text file, using a given precision
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <param name="iPrecBeforeDec"></param>
    ''' <param name="iPrecAfterDec"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                          ByVal sM() As Double, _
                          ByVal iPrecBeforeDec As Integer, _
                          ByVal iPrecAfterDec As Integer)
      'Save a matrix to file.
      Dim strF As String = ""
      Dim i As Integer

      Debug.WriteLine("Save numeric data to file " & strFile)
      If iPrecAfterDec = -1 Then
        strF = "0."
      Else
        For i = 1 To iPrecAfterDec
          strF = strF & "0"
        Next
        strF = strF & "."
      End If

      For i = 1 To iPrecBeforeDec
        strF = strF & "#"
      Next

      If System.IO.File.Exists(strFile) Then System.IO.File.Delete(strFile)

      Try
        FileOpen(5, strFile, OpenMode.Output, OpenAccess.Write)
        For i = 0 To sM.GetUpperBound(0)
          PrintLine(5, SPC(1), Format(sM(i), strF))
        Next i
        FileClose(5)
      Catch e As Exception
        MsgBox("Algebra.Save (file = " & strFile & "):" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
      End Try
    End Sub

    ''' <summary>
    ''' Save an integer vector to a file
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                                     ByVal sM() As Integer)
      'Save a matrix to file.
      Dim i As Integer

      Debug.WriteLine("Save numeric data to file " & strFile)
      If System.IO.File.Exists(strFile) Then System.IO.File.Delete(strFile)

      Try
        FileOpen(5, strFile, OpenMode.Output, OpenAccess.Write)
        For i = 0 To sM.GetUpperBound(0)
          PrintLine(5, SPC(1), sM(i).ToString)
        Next i
        FileClose(5)
      Catch e As Exception
        MsgBox("Algebra.Save (file = " & strFile & "):" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
      End Try
    End Sub

    ''' <summary>
    ''' Save an integer matrix to a file
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                                     ByVal sM(,) As Integer)
      'Save a matrix to file.
      Dim i, j, iHiRow, iHiCol As Integer

      If System.IO.File.Exists(strFile) Then System.IO.File.Delete(strFile)
      Try
        FileOpen(5, strFile, OpenMode.Output, OpenAccess.Write)
        iHiRow = sM.GetUpperBound(0)
        iHiCol = sM.GetUpperBound(1)
        For i = 0 To iHiRow
          For j = 0 To iHiCol - 1
            Print(5, sM(i, j).ToString(), SPC(1))
          Next j
          PrintLine(5, SPC(1), sM(i, iHiCol).ToString())
        Next i
        FileClose(5)
      Catch e As Exception
        MsgBox("Algebra.Save (file = " & strFile & "):" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
      End Try
    End Sub

    ''' <summary>
    ''' Save a byte vector to a file
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                                     ByVal sM() As Byte)
      'Save a matrix to file.
      Dim i As Integer
      Debug.WriteLine("Save numeric data to file " & strFile)
      If System.IO.File.Exists(strFile) Then System.IO.File.Delete(strFile)

      Try
        FileOpen(5, strFile, OpenMode.Output, OpenAccess.Write)
        For i = 0 To sM.GetUpperBound(0)
          PrintLine(5, SPC(1), sM(i).ToString)
        Next i
        FileClose(5)
      Catch e As Exception
        MsgBox("Algebra.Save (file = " & strFile & "):" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
      End Try
    End Sub

    ''' <summary>
    ''' Save a byte matrix to a file
    ''' </summary>
    ''' <param name="strFile"></param>
    ''' <param name="sM"></param>
    ''' <remarks></remarks>
    Public Overloads Shared Sub Save(ByVal strFile As String, _
                                     ByVal sM(,) As Byte)
      'Save a matrix to file.
      Dim i, j, iHiRow, iHiCol As Integer

      If System.IO.File.Exists(strFile) Then System.IO.File.Delete(strFile)
      Try
        FileOpen(5, strFile, OpenMode.Output, OpenAccess.Write)
        iHiRow = sM.GetUpperBound(0)
        iHiCol = sM.GetUpperBound(1)
        For i = 0 To iHiRow
          For j = 0 To iHiCol - 1
            Print(5, sM(i, j).ToString(), SPC(1))
          Next j
          PrintLine(5, SPC(1), sM(i, iHiCol).ToString())
        Next i
        FileClose(5)
      Catch e As Exception
        MsgBox("Algebra.Save (file = " & strFile & "):" & e.Message, CType(MsgBoxStyle.Critical + MsgBoxStyle.OkOnly, MsgBoxStyle))
      End Try
    End Sub

    Private Shared Sub GetBounds(ByVal sM(,) As Double, _
                                 ByRef iHiRow As Integer, _
                                 ByRef iHiCol As Integer)
      iHiRow = sM.GetUpperBound(0)
      iHiCol = sM.GetUpperBound(1)
    End Sub

    ''' <summary>
    ''' Create a string from a matrix
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function ToString(ByVal sM(,) As Double) As String
      Return ToString(sM, 99)
    End Function

    ''' <summary>
    ''' Create a string from a matrix, using a given precision
    ''' </summary>
    ''' <param name="sM"></param>
    ''' <param name="iNrAfterComma"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function ToString(ByVal sM(,) As Double, ByVal iNrAfterComma As Integer) As String
      Dim strText As String = vbNewLine
      Dim i, j, iHiRow, iHiCol As Integer
      Dim strFormat As String = "G" & iNrAfterComma

      GetBounds(sM, iHiRow, iHiCol)
      Dim strValues(iHiRow, iHiCol) As String, iMaxSize As Integer = 0

      For i = 0 To iHiRow
        For j = 0 To iHiCol
          strValues(i, j) = sM(i, j).ToString(strFormat)
          If iMaxSize < strValues(i, j).Length Then iMaxSize = strValues(i, j).Length
        Next j
      Next i

      iMaxSize += 1
      For i = 0 To iHiRow
        For j = 0 To iHiCol - 1
          strText = strText & strValues(i, j).PadLeft(iMaxSize)
        Next j
        strText = strText & strValues(i, iHiCol).PadLeft(iMaxSize) & vbNewLine
      Next i

      Return strText
    End Function

    ''' <summary>
    ''' Create a string from a vector
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function ToString(ByVal sV() As Double) As String
      Return ToString(sV, " ", 99)
    End Function

    ''' <summary>
    '''  Create a string from a vector, using a given precision
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="separator"></param>
    ''' <param name="iNrAfterComma"></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function ToString(ByVal sV() As Double, ByVal separator As String, ByVal iNrAfterComma As Integer) As String
      Return ToString(sV, separator, "G" & iNrAfterComma)
    End Function

    ''' <summary>
    '''  Create a string from a vector, using a given format string
    ''' </summary>
    ''' <param name="sV"></param>
    ''' <param name="separator"></param>
    ''' <param name="format "></param>
    ''' <returns></returns>
    ''' <remarks></remarks>
    Public Overloads Shared Function ToString(ByVal sV() As Double, ByVal separator As String, ByVal format As String) As String
      Dim strText As String = ""
      Dim i, iHiCol As Integer

      iHiCol = sV.GetUpperBound(0)
      For i = 0 To iHiCol - 1
        strText = strText & sV(i).ToString(format) & separator
      Next i
      strText = strText & sV(iHiCol).ToString(format)
      Return strText
    End Function

  End Class
End Namespace

